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ABSTRACT 

The non-linear regime of Ray leigh- Taylor instability (RTI) in a radiation supported atmosphere, 
consisting of two uniform fluids with different densities, is studied numerically. We perform simula- 
tions using our recently developed numerical algorithm for multi-dimensional radiation hydrodynamics 
based on a variable Eddington tensor as implemented in Athena, focusing on the regime where scatter- 
ing opacity greatly exceeds absorption opacity. We find that the radiation field can reduce the growth 
and mixing rate of RTI, but this reduction is only significant when radiation pressure significantly 
exceeds gas pressure. Small scale structures are also suppressed in this case. In the non-linear regime, 
dense fingers sink faster than rarefied bubbles can rise, leading to asymmetric structures about the 
interface. By comparing the calculations that use a variable Eddington tensor (VET) versus the Ed- 
dington approximation, we demonstrate that anisotropy in the radiation field can affect the non-linear 
development of RTI significantly. We also examine the disruption of a shell of cold gas being acceler- 
ated by strong radiation pressure, motivated by models of radiation driven outfiows in ultraluminous 
infrared galaxies. We find that when the growth rate of RTI is smaller than acceleration time scale, 
the amount of gas that would be pushed away by the radiation field is reduced due to RTI. 
Subject headings: methods: numerical — hydrodynamics — instability — radiative transfer 



1. INTRODUCTION 

The Rayleigh- Taylor Instability (RTI) occurs when 
a high density fluid is accelerated by a low den- 
sity fluid in a gravitational field or n et acceleration 
(e.g. IChandrasekhari [l96ll: IShard [1981 . In an astro- 
physical context, RTI with or without magnetic fleld 
has been stud i ed extensively both numerically (e.g. 
Jun et all [Tool IDimonte et all 120041: IStone fc Gardmerl 
2007allbri and exper imentally (e.g. lAiidrews fc Spalding 
ToonTTDalzieil [1991 . In most previous studies, sup- 
port against gravity was provided by a gradient of the 
gas pressure or magnetic pressure alone. In this pa- 
per, we are interested in the dynamics of the RTI when 
support against gravity is provided primarily by radi- 
ation pressure. This situation is expected, for exam- 
ple, in the inner region of accretion disks around com- 
pact object; at the interface of an H II region produced 
by massive star clusters and the surrounding medium 
(e.g. [Krumholz fc Matznei) l2009t iJacquet fc Krumholzl 
120111 ) and the dusty tori around lumi nous Active Galac- 
tic N uclei (AGN) (e.g. [Krolik 2007; ISchartmann et al.l 
I2010f ). In ultraluminous infrared galaxies (ULIRGS), 
radiation pressure on the dust grains may provide the 
dominant vertical support against gravity and acceler- 
ate the cold (neutral) gas outflows observed in these 
systems fe.g. lMurrav et all 120051: [Thompson et al.ll2005l: 
IMartinI I2005D . The structure of radiation dominated 
bubbles around super-Eddington massive stars also de- 
pend s on whether RTI dev e lops on the bubble interface 
(e.g. iKrumholz et~aL[ [2001 iKuiper et all [2012[ ) . Radi- 
ation RTI, together with Kelvin-Helmholtz instability, 
may also be important to generate the turbulence in 
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the i nterstellar medium (ISM) (e.g. IMcKee fc Ostrikeil 
l2007f ). Thus, understanding the effects of radiation fleld 
on linear growth rate and non-linear structures of RTI, 
as well as assessing the role of radiative RTI in driv- 
ing turbulence and possibly limiti ng the effectiveness 
of radiative support and dr iving (jKuiper et al.l 120111 : 
IKrumholz fc Thompsoi]l2012f ). is an important step to 
understand those physical systems. 

Effects of a radiation field on RTI have been 
studied analytically in the linear regime with dif- 
ferent approximations. In the optically thin limit, 
radiation flux only provi des a backgro u nd ac celer- 
ation, and both jMathcw s fc Blumenthal (|1977| ) and 
IJacquet fc Krumholzl (f20ll] ) conclude that stability de- 
pends only o n the direction of the effective gravitational 
acceleration. iKrolikI (|1977D studied the global RTI for 
an incompressible slab of gas accelerated by radiation 
pressure. Unstable mode exists if local radiative acceler- 
ation c orrelates positively w" i th tot al optical depth. Re- 
cently, IJacquet fc Krumhold (|201lD examined the effects 
of radiation field on the stability of the interface between 
two fluids in the isothermal and adiabatic limits. In the 
isothermal limit, the radiation flux only provides an ef- 
fective background acceleration. In the adiabatic limit, 
the radiation pressure plays the same role as gas pressure. 
However, the background state they envisage for the lin- 
ear analysis (two half infinite constant density planes) 
cannot precisely be in thermal equilibrium (see section 
[3|) with non-zero absorption opacity. This may affect the 
growth rate they calculate and the evaluation in the non- 
linear regime when the timescale for thermal evolution is 
comparable to or shorter than the growth rate of the 
RTI. 

If there is no density jump in a radiation supported at- 
mosphere without magnetic field, the nature of the sta- 
bility criterion is not completely settled. With the Ed- 
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dington approxiniation. [Spiegel &: Tad (|1999[ ) and lShavTvl 
(|2001f ) find that there are global hydrodynamic unstable 
modes for a hot atmosphere, even for luminosities moder- 
ately below t he Eddington fimit. However, the instabili- 
ties fou nd bvlSpieeel fc Tad (1999) a re not pr esent in the 
iShavivl (120011) analysis. Also. Spie gel fc Taol (l999) note 
that " iMarzed (flOTSfi arrived at a different conclusion by 
solving the full tran s fer eq uation numerically in a PhD 
thesis. [Turner et al.l ()2005l ) found overturning modes in 
simulations with a fixed temperature at the lower bound- 
ary, which are consistent with Shaviv's Type I modes, 
althoug h the author s find no evidence for the Type II 
modes (|Shavivll200il) . 

Here we will focus on the traditional case for RTI with 
a high density fluid on top of a low density fluid separated 
by an inflnite thin interface. In the previous analytical 
studies, simpliflcations are necessary to make the prob- 
lem tractable. For example, radiation pressure is usu- 
ally assumed to be isotropic, which is not true in general 
at the interface. In this work, we relax these assump- 
tions by solving the radiation hydrodynamic equations 
numerically. We study both the linear regime, using 
an Eddington tensor computed self-consistently from the 
time- independent radiation transfer equation (i.e., we al- 
low for anisotropic radiation pressure at the interface), 
and we also follow the RTI into the non-linear regime 
which is not possible in analytical studies. We note in 
passing that we have also used our numerical methods 
to test the stability of a radiation supported atmosphere 
with a smooth density profile , and find no evidence for 
instability, in agreement with [Turner et all (|2005). 

The paper is organized as follows. In section [21 we de- 
scribe the equations we solve and the numerical code we 
use. We then consider the problem of RTI with a single, 
initially static interface between two fluids of different 
density. We describe the background equilibrium state 
and initial perturbations in section[31 and summarize our 
results in section [Jj In section [SJ we describe simulations 
of RTI in shells being accelerated by radiation forces. We 
summarize and conclude in section [6] 



2. EQUATIONS 

We solve the radiation hydrodynamic equations in 
the mixed fr a me w ith radiation source terms given by 
[Lowrie et al.l (|1999D . We assume local thermal equi- 
librium (LTE) and that the Planck and energy mean 
absorption opacities are the same. Detaile d discussion 
of the equations we solve can be found in [Jiang et al.l 
(|2012f ). With a vertical gravitational acceleration g, the 
equations are 



dt 
dE 



| + V.(p.) = 0, 

V ■{pvv + P) = pg-^,{P), 



— + V-[{E + P)v]^pv ■ g - cSriE), 

-g^+V-Fr=cSr{E), 



Here, p is density, P = P\ (with I the unit tensor) and P is 
gas pressure, Ua and ag are the absorption and scattering 
opacities. Total opacity (attenuation coefficient) is at = 
o's + o'a- The total gas energy density is 



E^Eg + Ipv^, 



(2) 



where Eg is the internal gas energy density. We adopt 
an equation of state for an ideal gas with adiabatic index 
7, thus Eg = P/(7 - 1) and T = P/i?idoaiP, where i?idoai 
is the ideal gas constant. The radiation pressure is 
related to the radiation energy density E,. by the closure 
relation 



fE, 



(3) 



where f is the variable Eddington tensor (VET). Radia- 
tion fiux is Fr while c is the speed of light. The grav- 
itational acceleration g is fixed to be a constant value 
along the —z axis. Sr{P) and Sr{E) are the radiation 
momentum and energy source terms. 

Following [Jiang et al.l ([2012[ ). we use a dimensionless 
set of equations and variables in the remainder of this 
work. We convert the above set of equations to the di- 
mensionless form by choosing fiducial units for velocity, 
temperature and pressure as oq, Tq and Pq respectively. 
Then units for radiation energy density Er and flux Fr 
are a^TQ and carTg. In other words, = 1 in our units. 
The dimensionless speed of light is C = c/flQ. The orig- 
inal dimensional equations can then be written to the 
following dimensionless form 



d{pv) 

at 

dE 
'dt 



| + V.(p.) = 0, 

■{pvv + P)=pg-¥Sr{P), 
V ■[{E + P)v] =pvg- FCSr{E), 



dEr 
~dt 
dFr 
"df 



CV ■ Fr = CSr{E), 
CV- P^ = CSr(P), 



(4) 



1 dFr 

'^~dr 



V-P, = Sr(P). 



(1) 



where the dimensionless source terms are, 
Sr(P) = -a. +a.liT^~Erl 

Sr{E)=a,{T^ - Er) + (a, - a,)^ . {Pr ^^^L±^(|3) 

The dimensionless number P = arTp/ Po is a measure of 
the ratio between radiation pressure and gas pressure in 
the flducial units. We prefer the dimensionless equations 
because the dimensionless numbers, such as C and P, can 
quantitatively indic ate the imp o rtanc e of the radiation 
field as discussed in [Jiang et al.[ (|2012l ). 

We solve these equations in a 2D x — z plane with the 
recently developed radiation transfer module in Athena 
([Jiang et al.[ [2012[ ). The continuity equation, gas mo- 
mentum equation and gas energy equation are solved 
with modified Godunov method, which couples the stiff 
radiation source terms to the calculations of the Riemann 
fiuxes. The radiation subsystem for Er and Fr are solved 
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with a first order implicit Backward Euler raethod. De- 
tails on the n umerical alg o rithrn . and tests of the code are 
described in iJiang et aLl (|2012D . The variable Edding- 
ton tensor is computed from angular quadratures of the 
specific intensity J^, which is calculated from the time- 
independent transfer equation 



dir 



Kt{S - Ir). 



(6) 



Details on how we calculate t he VET, including tests, 
are given m IDavis et all (|2012D . 

3. BACKGROUND EQUILIBRIUM STATE 

As is usual, the background equilibrium state used to 
study RTI in this work is an interface which separates 
two uniform fluids with different densities. If there is 
no radiation field, gravitational acceleration (or effective 
gravitational acceleration) must be balanced by a gas 
pressure gradient. With radiation, the background state 
needs to satisfy both mechanical and thermal equilib- 
rium if the absorption opacity is nonzero. Because the 
background state is uniformly parallel to the interface, 
we only need to consider equation (j4]) perpendicular to 
the interface (and parallel to the direction of gravity), 
which is the z direction. The equilibrium state should 
satisfy the following equations 

dP 

— + pg = FatFr, 
oz 



dFr 
dz 



rp4 

= 0, 



(7) 



The second equation states that radiation temperature 
must be the same as gas temperature in thermal equilib- 
rium when there is a non-zero absorption opacity. The 
third equation means radiation flux is a constant. Comb- 
ing the four equations, and making the Eddington ap- 
proximation, we find 



-pg- 



(8) 



That is the gradient of the sum of gas and radiation 
pressure that balances gravity. The gas pressure P is 
related to the density p and temperature T via ideal gas 
equation of state 



P = RidcaipT = pT. 



(9) 



At the interface, if there is a density jump, there must 
be a corresponding jump in the gas temperature to sat- 
isfy continuity of the total pressure. In thermal equi- 
librium, the radiation temperature and gas temperature 
must be the same, which means that radiation temper- 
ature will also jump at the interface. In this case, the 
diffusion equation cannot be satisfied at the interface. 
That is, with radiation, and if absorption opacity is non- 
zero, this configuration cannot satisfy both mechanical 
and thermal equilibrium because gas pressure and radia- 
tion pressure depends on temperature in different ways. 
Thus an interface with a constant non-zero absorption 



opacity cannot be used as an equilibrium background 
state. 

One way of circumventing this constraint is to sim- 
ply consider a background that is not strictly in thermal 
equilibrium, but with a thermal timescale that is much 
longer than the growth rate of the instabilities of inter- 
est. One can then assume that the evolution of the back- 
ground has only minor effects on the stability properties 
and subsequent evolution. The energy equation in (j4]) 
indicates that the timescale to achieve thermal equilib- 
rium can be reduced by making the flow less relativistic 
(lower C), less radiation dominated (lower P), or by low- 
ering the absorption opacity. If we take the limit that 
absorption opacity goes to zero on one (or both) sides of 
the interface, and T decouple and T can be discontin- 
uous while Er remains continuous. From this point on, 
we consider domains where absorption opacity is zero 
everywhere, and the thermal time is effectively inflnite. 

3.1. Equilibrium State with Pure Scattering Opacity 

In order to study the evolution of density discontinu- 
ities due to RTI for radiation pressure supported inter- 
face, and to compare to previous studies of RTI of inter- 
faces in gas pressure supported atmosphere, we only con- 
sider material with pure scattering opacity. The specific 
scattering opacity k is assumed to be a constant (as in the 
case of electron scattering opacity). Then <Jt = pK, and 
radiation temperature can be independent of gas tem- 
perature. 

The background equilibrium state with pure scattering 
opacity is calculated in the following way. For a given 
density, opacity and gravitational acceleration, many 
equilibrium states can be constructed, depends on the 
relative contributions from the gas pressure and radia- 
tion pressure to support the gravity. We first choose a 
constant flux Fr o and deflne 



r,0 



pg 



(10) 



This a is the fraction of gravitation acceleration that is 
balanced by radiation pressure gradient and 1 — a is the 
fraction balanced by gas pressure gradient. The radiation 
pressure proflle can then be calculated from the diffusion 
equation, and the gas pressure proflle can be calculated 
according to dP/dz = — (1 — a)pg. Solutions in the two 
fluids are matched in the interface by requiring that gas 
and radiation pressure are continuous. 

There are two special background equilibrium states. 
The flrst case is a = 0, which means that the initial back- 
ground radiation flux Fr^o is zero and gravity is balanced 
by gas pressure gradient alone. This is very similar to 
previous studies of RTI except that now the material is 
embedded in a uniform radiation held. The second case 
is a = 1, which means that the gas pressure is uniform 
and gravitational acceleration is balanced by radiation 
pressure gradient alone (Eddington limit). We will focus 
on the two special cases flrst. 

3.2. Simulation Setup 

The 2D simulations are done in the (x, z) plane with 
gravitational acceleration along the —z direction. In our 
units, the size of the simulation box is (—0.5, 0.5) x (—1, 1) 
and a resolution of 128 x 512 grid points is used. The 
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Table 1 

Parameters of the simulations with only one interface. 
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Rad 
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Single 
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Rad 
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Rad 


1 
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16 


1 


Single 
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Hydro 
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Random 


G 


Rad 
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Random 
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Rad 
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Rad 
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1 
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eigenvector 
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Rad 


10^ 


10^ 
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eigenvector 
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Rad 
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eigenvector 



gravitational acceleration is g — 0.1, the dimensionless 
speed of light C is 10^, while gas pressure P = 10 and 
Er — 10 at the interface z = 0, except for simulations L 
— P, where we use P = 1 and Er = 3 a.t the interface. 
The density in the upper region z > is denoted by p+ 
while it is p_ in the bottom half z < 0. The isothermal 
sound speed in p_ is Cg = 1. Time is measured in units 
of the sound crossing time along the horizontal direction. 
The ratio between the free-fall velocity to the isothermal 
sound speed is {gLz)^-^/cs — 0.45. 

Periodic boundary conditions are used in the horizon- 
tal direction, while the gradient of both gas and radi- 
ation pressure in the vertical direction are extended to 
the ghost zones to achieve a better hydrostatic equilib- 
rium. Radiation flux along the z direction is fixed to 
be the initial value to balance the gravity in the ghost 
zones while the x component is continuously copied to 
the ghost zones. The boundary condition of the short 
characteristic module used to calculate the VET is set 
such that a constant flux is also maintained in that mod- 
ule. Since the divergence of the radiative flux is zero due 
to the assumption of negligible absorption opacity, this 
simply means we adjust the incoming radiation field at 
the base of the domain to provide the appropriate flux at 
the bottom boundary. For simplicity we assume isotropic 
radiation for the incoming intensity, although the angular 
distribution will generally be problem dependent. This 
allows us to ensure that the ratio of radiation flux to radi- 
ation energy density from the radiation transfer module 
will be consistent with the result obtained by evolving 
eqs. (|3]). If this consistency is not guaranteed, the incor- 
rect VET can cause numerical instability near a stable 
interface. 

The underlying issue is that the boundary conditions 
on the short characteristics solver should be consistent 
with the initial condition and boundary condition on the 
radiation energy and momentum equations. For exam- 
ple, if the domain is assumed to be embedded in an opti- 
cally thick background so that <C Er initially, then the 
Eddington tensor should be nearly uniform (~ 1/3) ev- 
erywhere. If the boundary conditions on the short char- 
acteristics method due not make the same assumption 
(e.g. if they assumed the upper boundary corresponded 



to vacuum), the VET would increase appreciably from 
1/3 near the upper boundary even though it should be 
constant. In this case, very small variations in Er can 
drive large and unphysical variations in Fr- We find that 
forcing the boundary conditions in the modules to return 
a consistent ratio of flux to energy density ensures a con- 
sistent Eddington tensor and ameliorates this problem. 

We consider three sets of perturbations to the initial 
background. In the first two sets, only the vertical ve- 
locity is perturbed, either by a single mode in the form 

= 0.0025(1 -I- cos(27r2;))(l -I- cos(7rz)), (11) 

or by a random collection of modes in the form 

= t'z,o(l + cos(7rz)), (12) 

where Vzfi is a random number uniformly distributed 
between —0.025 and 0.025. We also consider another 
set of single mode perturbations where we perturb 
Vx,Vz, Pr, Frx, and Frz based on the linear analysis de- 
scribed in the appendix. We normalize Vz as 

{vz)± = 0.01 cos(27rx) exp(=F27rz), (13) 

with perturbations in other quantities chosen to corre- 
spond to the incompressible eigenvectors. Other param- 
eters of the simulations are listed in Table [1] We only 
present results for interfaces with > p_. As in the 
purely hydrodynamic case, the interface is expected to 
be stable when the lower density fluid overlies the higher 
density fluid. We have confirmed this numerically for 
several test cases. 

4. RESULTS 

In this section, we describe the linear growth and non- 
linear evolution of RTI for different background states. 
By comparing a set of controlled simulations, we explore 
the behavior of RTI with different parameters, such as 
opacity and the ratio between radiation pressure and gas 
pressure for the two background states a = and a — 
1. We will first use the Eddington approximation by 
assuming radiation pressure to be isotropic, thereafter 
f = 1/31. We then relax this assumption and study the 
effects of anisotropic radiation pressure on RTI. 

4.1. Growth rate for a single mode perturbation 

For a non-radiative, gas pressure supported inter- 
face, the linear regim e has been studied extensively (e.g. 
IChandrasekhail [19611 ) . For a mode with wavelength A 
and thus wavenumber k = 2tt/X, the growth rate n of 
the inviscid, incompressible RTI is given by 

n=.Lk^±—P-. (14) 

However, for the radiation pressure supported interface, 
it is not possible to derive a general, but still simple for- 
mula for the growth rate of the radiation RTI. Nonethe- 
less, we can still examine the behavior of RTI in vari - 
ous regimes. For example, IJacauet fc Krumhold (|2011[ ) 
derived the growth rate of RTI for the radiation sup- 
ported interface in the optically thin as well as adiabatic 
limit. However, it is likely that their results in the adia- 
batic limit will be affected by their non-equilibrium back- 
ground state. In the optically thin limit, the background 
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radiation flux plays the role as an effective gravitational 
acceleration and so the growth rate will be very similar to 
the value given by equation ([M]) with g replaced by the 
effective gravitational acceleration g^s- Since t/eff = for 
a = 1, this would correspond to a marginally stable state 
with zero growth rate. In the very optically thick limit 
when photons are fully coupled with the gas and radia- 
tion field is isotropic, the radiation pressure behaves in 
the same way as gas pressure and so the growth rate will 
be very similar as given by equation (|14|). (See the ap- 
pendix for further discussion.) The intermediate regime 
is the most complicated case, as the Eddington tensor 
varies with z and standard linear analysis if of limited 
utility. 
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Figure 1. Left: Growth of the vertical displacement Zmax between 
the top of the low density fluid bubbles and the initial interface due 
to RTI for different bac kground states with a single mode pertur- 
bation given by eq. Illip . Simulation parameters are listed in Table 
[T] The black and blue short dashed lines are the theoretical pre- 
dictions from eq. I I14I I for interfaces with density ratios / p_ of 
4 and 16, respectively. Right: The same as the left panel, except 
for simulations initially perturbed with the eigenvectors. 

In addition to the density ratio at the interface, we an- 
ticipate that the behavior of RTI for radiation supported 
interface could differ depending on the background ra- 
diation properties, such as the ratio between radiation 
pressure and gas pressure and the characteristic optical 
depth. We proceed to vary these parameters and measure 
the growth rate from our simulations. We quantify the 
growth rate of RTI using 2:niaxi the maximum distance 
of the interface from z — 0. Zmax is measured by the 
perturbed regions that have density different from the 
initial density in those positions by at least 10%, such as 
the bubbles and fingers. 

We first consider the case with single mode eigenvector 
perturbations. As discussed in the Appendix, we expect 
the growth rate consistent with eq. ([T^ since C ^ 1 
and the Eddington tensor is diagonal. The results are 
shown in the right panel of Figure [T] For a density ratio 
p+/p- =4 and wavenumber k = 2n, eq. ([T4|) predicts 
■Zmax oc exp(0.61t), which is shown as solid line in the 
right panel of Figure [1] We find initial growth is faster 
than predicted for Zmax ^ 0.02 (corresponding to < 5 
cells). At such early times the vertical deformation of 
the interface is not well resolved and increases in the 
resolution improve agreement with constant growth rate 
solution. At later times, non-linear effects become im- 



portant and the growth rate drops below the prediction 
for all runs, including the purely hydrodynamic case. 

Although the curves are not uniformly consistent with 
exponential growth at a single constant growth rate, the 
range of growth rates inferred from the right panel of 
Figure [1] are close to the prediction of linear theory. Fur- 
thermore, many aspects of the linear theory are borne out 
by the results. Since P and k do not appear in eq. (fMl) 
the linear theory predicts no dependence on these param- 
eters to first order in C~^. This is essentially correct for 
the initial stages of growth where the solutions remain 
in the linear regime. However, after about one sound 
crossing time in the simulations with large P and k (iV, 
O, and P), the growth rate of Zmax slows relative to the 
hydrodynamic and low P cases. This appears to be a 
non-linear behavior as the amplitude of the perturbation 
grows. 

We attribute this non-linear behavior primarily to the 
effects of radiation drag, by which we mean that the 
fluid motions are damped by oppositely directed radi- 
ation forces. This requires that velocity fluctuations (on 
average) tend to directed opposite to the comoving radi- 
ation flux. Note that this is not equivalent to radiation 
viscosity, which is absent due to our assumption tha t the 
radiation field is isotropic (|Mihalas fc Mihalaslll984f) . To 
linear order in perturbed quantities and lowest order in 
C~^, radiation drag is only relevant for compressible mo- 
tions. The effects of radiation drag are most apparent in 
simulation where both P and k arc large. Runs O and 
P, with different values P and k but the same product of 
Pk, follow similar trajectories. This is consistent with the 
linear analysis since the characteristic frequency v oc Pk 
controls the importance of radiation drag on the solution. 
See the Appendix for further discussion. 

The results of right panel of Figure [1] can be compared 
with the second set of single mode perturbations that 
are not eigenvectors, which are also shown in the left 
panel of Figure [TJ In this similar diversity of 

growth rates is observed, with notably reduced growth 
for the runs with larger P and k. As with the eigenvec- 
tor perturbations this trend of slower growth with larger 
radiation pressure appears to be a result of the radia- 
tion drag. In this case, the initial perturbations give rise 
to compressive motions which are largely absent in the 
eigenvector initial conditions. This allows radiation drag 
effects to become important at lower amplitude than in 
the eigenvector runs. But the range of late time (non- 
linear) evolution observed is qualitatively similar in the 
two sets of runs. 

We have performed another simulation with the same 
parameters as simulation D, but with a = i.e. with a 
background state supported by a gas pressure gradient 
instead of a radiation pressure gradient. We find a very 
similar growth history to that in simulation D. This con- 
firms that the difference in the role of gas and radiation 
pressure in supporting the background has no qualitative 
effect on the resulting behavior. Radiation has relatively 
modest effects when P is small, but generically leads to 
slower growth when P and n are large. This strengthens 
the conclusion that damping of the compressible motions 
by radiation drag is the main effect of the radiation field. 
The fluid simply has more difficulty moving with respect 
to strong radiation field. 
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4.2. Non-linear structure with a single mode 
perturbation 

The non-linear structure developed by RTI with single 
mode initial perturbations are shown in Figure [H taken 
at late times in simulations A — D. By comparing the left 
two panels (simulation A and C), we see that when the 
radiation pressure is comparable to the gas pressure and 
radiation pressure is assumed to be isotropic, RTI for the 
radiation supported interface and gas pressure supported 
interface have very similar non-linear structures. This is 
consistent with the fact that they have very similar linear 
growth phase as shown in Figure [U especially when 
is increased to 16 (simulation E). 

The non-linear structures for the large radiation pres- 
sure cases (simulation B and D) are quite different from 
simulations A and C, especially at small scales. The 
Kelvin-Helmholtz (KH) instability, which is produced 
during the non-linear phase of RTI because of the shear- 
ing between the rising bubbles and falling fingers, is 
clearly damped in simulation B where the optical depth 
across the box is 1. When optical depth is increased to 
10^ for simulation D, the small scale structures can sur- 
vive because the size of those structures are much larger 
than photon diffusion length. However, these structures 
are still quite different when compared to the struc- 
tures in simulations A and C. The vortices are strongly 
stretched along the vertically direction when a strong ra- 
diation pressure exists. 

4.3. Non-linear structures of RTI with random initial 
perturbations 

To study the interactions between different modes of 
RTI, we add random perturbations to the vertical veloc- 
ity as given by equation ([T2|) . As in section 221 the effect 
of radiation field on the RTI can be shown by comparing 
a set of controlled simulations {F, G, H, /, J), the pa- 
rameters of which are listed in Table [TJ The Eddington 
approximation is still used in those simulations. 

In Figure |3l we show four snapshots from simulations 
F, G, H, J (from left to right and from top to bottom) to 
illustrate the range of non-linear structures that we find. 
By comparing the top two panels {F and G) , we see that 
when radiation pressure is comparable to gas pressure 
(P = 1), the non- linear outcome of RTI for the radia- 
tion supported interface is similar to the case with gas 
pressure supported interface. There are two major differ- 
ences when a radiation field is present. The first is that 
the dense, heavy fingers sink faster than the light bubbles 
rise up. This is not the case for the non-radiative RTI. 
Another difference is that for radiative RTI, the fluid is 
less mixed. In non-radiative RTI, there are many cells 
with density distributed between the whole range of the 
maximum and minimum density due to the small scale 
turbulence driven by secondary KH instabiity. However, 
with radiation field, there are fewer such cells. The re- 
sult when /9+ is increased by a factor of 4 is shown in 
the bottom right panel (simulation J). The mixing rate 
is increased as cells with intermediate densities occupy a 
larger fraction of the volume. The non-linear structures 
in this case bear more similarity to those in the purely 
hydrodynamic run. Recall that simulation E has a lin- 
ear growth rate very close to the growth rate of normal 
RTI, as shown in Figure [TJ These results suggest that 




Figure 2. Density and velocity field for the non-linear regime of 
RTI with a single mode perturbation. From left to right, top to 
bottom, the simulations are A at time 8.8, B at time 14.0, C at 
time 11.6 and D at time 29.3. Parameters of the simulations are 
listed in Table[T] The growth history of the four cases are shown in 
Figure [T] All the simulations adopt the Eddington approximation 
f = 1/31, which means the radiation pressure is assumed to be 
isotropic. 




Figure 3. Density and velocity field for the non-linear regime of 
RTI with random initial perturbations. From left to right, top to 
bottom, the simulations are F at time 15.2, G at time 15.8. H at 
time 26.0. J at time 10.7. All the simulations use the Eddington 
approximation. Parameters of the simulations are listed in Table 

m 
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a higher ratio of pj^j p- reduces the impact of radiation 
and leads to behavior more consistent with the purely 
hydrodynamic RTI. 

The effects of strong radiation pressure (P = 10"*, sim- 
ulation H) on the non-linear structure is shown in the 
left bottom panel of Figure [31 where optical depth across 
the box along the horizontal direction is 1. Due to the 
strong damping effect of the radiation field, all the small 
scale perturbations are damped and only the mode with 
the longest horizontal wavelength grows. In other words, 
the radiation field sets a minimum scale, above which the 
turbulence due to RTI can be sustained when radiation 
pressure is significant. Then the non-linear outcome is 
very similar to that in the runs with single mode pertur- 
bation with the same wavelength (simulation B shown in 
Figure [5]) . The bubbles and fingers are well defined and 
mixing ratio is significantly reduced. 

As in section 14.11 we use ^max and Zmin to quantify 
the properties of RTI for the case with initial random 
perturbations. Here Zmin is similar to Zmax but only for 
the bottom half space (z < 0) while Zmax is only for the 
top half space {z > 0). Growth history for the five sim- 
ulations {F, G, H, I, J) are shown in Figure SI For non- 
radiative RTI (simulation F) , the structure is symmetric 
with respect to the initial interface, which means Zmax 
and Zmin grow at approximately the same rate. However, 
this is not true for the RTI when a strong radiation field 
is present. As shown in Figure IH for the case with radi- 
ation supported interface, Zmin reaches —1 at an earlier 
time than the time when Zmax reaches 1 . The asymmetry 
is reduced when opacity is increased (simulation /) . The 
non-monotonic growth of Zmax and Zmin is likely caused 
by the mixing of the cells at the top of the bubbles and 
fingers. When the mixing is reduced (as in simulation 
H), the changes of Zmax and Zaiin are always monotonic. 

In the left panel of Figure [SJ we show the net accel- 
eration due to gravity and radiation fiux along the ver- 
tical direction, Onet = PK-FrCz — 5, for simulation G at 
time 15.8. The initial hydrostatic equilibrium state has 
Onet = everywhere. The snapshot shows that for the 
bubbles moving up, a,ict > 0, which means that they are 
pushed up by the radiation force. For the falling fingers, 
flnct < 0, which means that radiation force is not able to 
balance the gravity. The right panel of Figure [5] shows 
the history of the change of radiation energy density dEr 
and the divergence of the radiation fiux across the whole 
simulation box dF,. with respect to the initial values. De- 
pending on whether the density in each cell is larger than 
0.5(p+ + P-) or not, we can locate the regions contain- 
ing high density fluid or low density fluid. The radiation 
energy density inside the high density fluid is increased 
because the high density fluid falls from the low radiation 
energy region to the high radiation energy region. The 
situation is inverted for the low density fluid. Because 
of the density (and, therefore, photon mean free path) 
difference, the radiation energy density of the low den- 
sity fluid is changed more quickly than the high density 
fluid. Thus the bubbles lose the support from radiation 
force more quickly than the fingers and the asymmetry 
arises. When opacity is increased for the fingers, the pho- 
ton mean free path is reduced and the diffusion time is 
longer. 

When density contrast is increased by a factor of 4 
(simulation J), the growth rate is increased and the 




Figure 4. Growth of vertical displacements Zmax and Zmin (the 
largest distance between the initial interface and the perturbed 
regions for 2; > and z < respectively) for RTI with random 
perturbations. The black (red) lines are for Zmax (^min)- Evolution 
over a longer time scale for simulation H is shown in the left corner. 
Parameters of the simulations are given in Table [l] 




Figure 5. Left: Snapshot of velocity (the vectors) and net accel- 
eration (the color) due to radiation and gravity at time t = 15.8 
for simulation G. The corresponding density distribution at this 
time is shown in the top right panel of Figure [S] Right: The top 
panel shows the temporal evolution of the change of radiation en- 
ergy density for low density fluid (dashed black line), high density 
fluid (solid black line) and total radiation energy density (red line) 
for simulation G. The bottom panel shows the change of the diver- 
gence of the flux across the whole simulation box compared with 
the initial value. 
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asymmetry is decreased. When radiation pressure is 
large enough (P = 10* for simulation H), the growths of 
2^max and Zmin are significantly decreased, which is consis- 
tent with what we find for the single mode perturbation 
shown in section |43] 

4.4. Mixing of the two fluids in the Rayleigh- Taylor 
Instability 

With random perturbations, the non-linear regime of 
the RTI generates turbulence, in which there is signifi- 
cant mixing. 

To quantify the degree of mixing, we count the num- 
ber of cells with density between 90% of the maximum 
density d+ and 1.1 times the minimum density d-. The 
degree of mixing can be estimated by the ratio between 
the mixed cells and the total number of cells. The mix- 
ing ratio for the simulations with initial random pertur- 
bations is shown in Figure [H For simulation G with 
K = 1 and comparable gas pressure and radiation pres- 
sure, the mixing ratio is only 71% of the mixing ratio 
for non-radiative RTI. When we increase the opacity to 
K = 10^ (simulation /), the mixing ratio is very close to 
the case of non-radiative RTI because the photons are 
so tightly coupled to the fluid. When all the available 
modes are optically thick, the mixing that results from 
RTI is very similar in the radiative and non-radiative 
cases. However, if the smallest eddies are optically thin 
while the largest modes are optically thick, photons will 
decouple from the fluid on small scales and there will 
be less mixing on these scales. When radiation pressure 
is increased to P = 10"*, the mixing ratio is dropped to 
~ 8%. This is consistent with the fact that all the small 
scale structures are damped and KH instability, which 
is the most important cause of mixing the two fluids in 
normal RTI, is suppressed by the strong radiation field. 
We caution that we have not performed a numerical con- 
vergence study of mixing in the radiation RTI, and the 
precise mixing fractions we quote in this section are likely 
to change with numerical resolution. In fact, without ex- 
plicit viscosity, converged results for the mixing rate are 
not possible. Nonetheless, the trends in the mixing rate 
with radiation pressure and opacity reported in this sec- 
tion should be independent of the viscosity and therefore 
numerical resolution. 

4.5. Effects of Anisotropic Radiation Pressure 

All the simulations listed in Table [1] adopt the Edding- 
ton approximation for the radiation field. Since the ra- 
diation pressure is assumed to be isotropic, these results 
only apply to very optically thick flows. We now consider 
cases with moderate to low optical depth and take into 
account the anisotropy of the radiation field by calculat- 
ing VET self-consistently with equation (|6]). This allows 
us to assess the impact of anisotropic radiation pressure 
on the development of the RTI. Here we focus on the 
cases with a constant background radiation flux. 

We first carry out a simulation G2 with Eddington 
approximation, which has almost the same setup as sim- 
ulation G with a different initial perturbation. Here we 
perturb the density randomly as 

p = Pa + 6p{l + cos^nz)), (15) 

where Sp is a random number uniformly distributed be- 
tween — 0.5po and 0.5po and the perturbation is only ap- 



plied when \z\ < 0.5. The initial radiation energy density 
at the bottom is 1.72 and decreases to 0.22 at the top due 
to the background radiation flux -FV.Oi which balances the 
gravity. All the other initial conditions are the same as 
simulation G. A third simulation employing the VET is 
labeled as G3. Incoming specific intensities are applied 
at both the top and bottom of the box to maintain the 
same constant background flux i^r,o- 

The initial vertical profile of the Eddington tensor is 
shown in Figure [T] For the pure scattering opacity case, 
the anisotropy is mainly caused by the stratification of 
the radiation energy density in order to provide the back- 
ground radiation flux. From bottom to top, f^x decreases 
while fzz increases monotonically. The Eddington ten- 
sor changes vary rapidly within roughly one optical depth 
from the boundary and is close to 1/3 in the middle of 
the simulation box. If we increase the vertical size of 
the simulation box, the profile of the Eddington tensor 
within one optical depth from the boundary remains the 
same while the middle region will become more isotropic. 
Not that we do not find fxx — fzz = 1/3 at the bottom 
boundary because the downward going (backscattered) 
radiation field is not perfectly isotropic. 

The growth history and mixing ratio of the two simu- 
lation G2 and G3 are shown in Figure [51 It is interesting 
to see that growth rate of G3 with VET is slower than 
the growth rate of G2. The mixing ratio from G3 is also 
smaller than the mixing ratio from G2. In Figure |9l we 
show snapshots of the non-linear structures to illustrate 
the differences between the two simulations. Consistent 
with Figure |H1 there are more mixed (intermediate den- 
sity) cells from simulation G2 and fewer in simulation G3. 
The top and middle panels in Figure |9] also show that the 
characteristic scale of eddies in G3 is larger than that of 
simulation G2. The slower growth rate in G3 is consis- 
tent with equation (|14|) . which predicts that the growth 
rate of RTI decreases with decreasing wave number. 

Components of the Eddington tensor fxx and fzz at 
the same time for simulation G3 are also shown in the 
two bottom panels of Figure [9l Interestingly, fxx and 
fzz actually have similar structures to the density dis- 
tribution, which is also the photon mean free path dis- 
tribution. There is a clear transition for the Eddington 
tensor between the fingers and bubbles: fxx is enhanced 
inside the finger while fzz is decreased compared with 
the initial values. Considering the complicated behavior 
of the Eddington tensor with position, it is not surprising 
that the non-linear structures produced by the RTI are 
different for simulations with VET and the Eddington 
approximation. Moreover, this result is a cause for con- 
cern regarding the reliability of results computed with 
much more approximate methods, such as flux-limited 
diffusion. 

5. RADIATION SUPPORTED COLD SHELL 

In ULIRGS, radiation from OB association is 
usuall y super-Eddin g ton for adjacent dus ty gas 
(e.g. 'Murrav et al.' '2005!; [Thompson et all 120051 : 
iKrumholz fc Matzner 2009f ). The radiation field 
from these stars may have important feedback effects 
on the interstellar mediu m, such as reducing the star 
formation efficiency (e.g. [Thompson et al.l I2005D and 
driving the cold neutral outflows seen in many of these 
systems (|Martin ,2005, ). Momentum driving may be 
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Figure 6. Mixing ratio for the simulations with initial random 
perturbations. Parameters of the simulations are given in Tabic [l] 
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Figure 7. Vertical profile of the initial x — x (fxx) and z — z 
ifzz) components o f the Eddington tensor for the simulation G3 
described in section [4.51 The Eddington tensor is invariant along 
the X axis due to the symmetry of the initial condition, fxx de- 
creases and fzz increases monotonically from the bottom to the 
top. 





Figure 8. Effects of VET on the growth rate and mixing ra- 
tio for RTI in a radiation supported interface. The left panel 
shows the evolution of Zmax and Zmin as in Figure |4] while the 
right panel shows the mixing ratio as in Figure Simulation G2 
adopts the Eddington approximation, while a VET is calculated 
self-consistently for simulation G3. All the other initial conditions 
are the same for the two simulations. 



more effective than thermal driving at accelerating 
clouds to escape velocities while not overheating and 
ablating the cold neutral gas (e. g. iMurrav et al.l 120051 : 
[Martin 200S iSocrates et al1[20^ 

It has been recently argued |Krumholz fc Thompson! 
I2012f l that the RTI limits the ability of radiation pressure 
to support ULIRGs against gravity due to the tempera- 
ture sensitivity of the dust opacity. Here we consider a 
related but somewhat different question of whether radi- 
ation pressure can effectively accelerate cold shells of gas 
before the RTI grows to a point that it limits radiative 
driving or provides enough mixing that the shell is no 
longer cold relative to the ambient medium. 

Here we give a simplified two dimensional model to 
demonstrate the effects of RTI on the gas in the inter- 
stellar medium. Initially, a dense cold shell is in gas 
pressure equilibrium with surrounding hot gas. Grav- 
ity is balanced by the radiation pressure gradient. We 
assume pure constant scattering opacity so that we can 
start from an equilibrium state. Note that real dust can 
both absorb and scatter photons, with absorption opacity 
comparable to or slightly larger than scattering opacity 
at the temperature relevant for ULIRGS. 

W e adopt the parameters characteristic for ULIRGS 
(e.g. [Thompson et al.|[2005HMurrav et al.ll2005D . subject 
to some constraints on the dynamical range that can be 
efficiently evolved by our simulation methods. The num- 
ber density of the cold medium is ric — 1000 cm~'^, which 
corresponds to a mass density pc ~ 3.3 x 10^^^ g cm~^. 
The density ratio between the cold and hot medium is 
100. Temperature of the cold medium is Tc = 50 K. The 
cold and hot gas are in gas pressure equilibrium. Grav- 
itational acceleration is assumed to be 17 = 10~® cm/s^. 
A typical size scale is chosen to be /q ~ 1 pc- With these 
parameters, the sound crossing time for Iq is ~ 10^ yr. To 
see the effects of varying the opacity, we try two differ- 
ent scattering opacities 1 cm^ g~^ and 100 cm~^ g~^ so 
that the optical depths across for the cold gas, which 
is labeled as t^, are 0.01 and 1 respectively. Background 
radiation fiux is chosen such that acceleration due to ra- 
diation force is 1.01 times the gravitational acceleration. 
Thus, the cold shell will be accelerated upwards, repre- 
senting the dusty gas which is pushed away by the strong 
radiation field. Taking pc,Tc and Iq as density, tempera- 
ture and length units in our simulation, the dimensionless 
parameters P = 8.1 x 10^ and C = 7.10 x 10^. 

The size of the simulation box is x 4^o and the cold 
shell is located at the middle, which extends from — 0.5Zo 
to 0.5Zo- Periodic boundary conditions are used along the 
horizontal direction. For the vertical directions, density 
is fixed to be O.Olpc in the ghost zones. The vertical 
components of the velocity and radiation fiux in the ghost 
zones are fixed to the initial values while the horizontal 
components are copied from the last active zones. This 
vertical boundary condition can maintain the constant 
radiation flux entering and leaving the simulation box. 
Numerical resolution is 128 x 512 for all the simulations. 

As discussed in section l475l to calculate the VET cor- 
rectly, we find it is necessary to make the ratios of radi- 
ation fiux to radiation energy density from the radiation 
transfer modules to be consistent with that obtained by 
evolving eqs. The boundary conditions on the ra- 
diative transfer model enforce a constant flux and as- 
sume that the incoming intensity at the lower boundary 



10 




Figure 9. Comparison of tlie non-linear structure from simula- 
tions G2 and G3, sliowing tlie effeet of anisotropy in the radiation 
field captured through use of a VET. Top: Density distribution 
for G2 at time 15.2 (left) and 21.6 (right). Middle: Density dis- 
tribution for G3 at time 21.6 (left) and 28.9 (right). The vectors 
in the those panels are the radiation flux with the background flux 
subtracted. Bottom: The horizontal (fxx) and vertical (fzz) com- 
ponents of the VET for simulation G3 at time 28.9. 
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Figure 10. Initial profiles of the x — x and z — z components of 
the Eddington tensor for the radiation supported cold shell shown 
in Sectionis] The red lines are for the case when the optical depth 
across the cold shell is 1 while the black lines are for the case with 
opacity 100 times smaller. The rapid change in the Eddington 
tensor at positions z = —0.5 and z = 0.5 are due to the jump in 
density (and, therefore, photon mean free path) at the interfaces. 
Although the density contrasts are 100 for both the red and black 
lines, the change of Eddington tensor is much larger for Tj, = 1. 



is isotropic. 

Given the density distribution and opacity, the initial 
VET is shown in Figure [TO] for the two cases Ts = 0.01 
and Ts = 1. When Ts = 1, the cold shell is optically thick 
while the hot ambient medium is optically thin. Thus 
we see a rapid change of Eddington tensor near the in- 
terfaces at z — ±0.5/o. The change is much smaller when 

= 0.01 and the radiation field remains nearly isotropic, 
consistent with the isotropy of the incoming radiation at 
the base of the domain. The vertical component of the 
Eddington tensor fzz increases with height while the hor- 
izontal component f^x decreases with height. Anisotropy 
of the radiation field is larger when Ts = 1 and much 
smaller when = 0.01. 

Initially, we perturb the density randomly according 
to equation [151 The random number Sp is uniformly dis- 
tributed between —0.005 and 0.005. As the acceleration 
due to the photons is 1% larger than the gravitational 
acceleration in this setup, the whole shell will move up- 
wards while RTI develops. The structure of the cold shell 
at times 2.3 x 10^ yr and 1.1 x lO'' yr for the case Ts — I 
are shown in Figure [TT] RTI develops quickly at the 
lower interface, where the high density medium is on top 
of the low density medium. The finger makes the cold 
dense gas fall back towards the low density gas exponen- 
tially, instead of being pushed away by the photons. In 
the top panel of Figure [121 we show how the fraction of 
cold gas that is below the line —0.5lo, /mass, change with 
time. Initially all the cold gas is located between ~Q.5la 
and 0.5^0 and the whole shell is being pushed upwards 
before RTI develops, therefore /mass remains zero when 
t < 3.6. After RTI grows significantly, /mass increases to 
13% after just 2.8 x 10^ yr. The maximum and minimum 
heights of the cold dense gas, Zmax and Zmin, are shown 
at the bottom panel of Figure [T^ Consistently with the 
history of /mass, the whole shell moves upwards together 
during the time t < 3.2. After t = 3.2, Zmin decreases 
dramatically instead of increasing with time. The cold 
gas eventually reaches the bottom of the simulation box. 

To see the effects of opacity, we decrease the opacity by 
a factor of 100 so that — 0.01. In Figure [T51 we show 
the density and velocity distribution at times 2.3 x 10^ yr 
and 1.1 X 10^ yr for this case. The evolution of the cold 
shell is quite different in this case compared with the case 
Ts = 1, in that the two interfaces are both stable during 
the time when the shell crosses the simulation box. The 
whole shell is just pushed away by the photons and no 
cold gas is left. 

By comparing Figure [TT] and Figure [T51 we can also see 
the different effect of radiation force in the optically thick 
and thin regimes. Initially in both cases, a constant radi- 
ation acceleration, which is larger than the gravitational 
acceleration, is applied across the shell to push it up- 
wards. In the optically thin regime, the acceleration due 
to the radiation force is almost unchanged and it behaves 
like an effective gravitational acceleration. However, in 
the optically thick regime when RTI develops, the radi- 
ation flux is enhanced inside the bubbles and reduced 
inside the fingers. In this case the whole shell no longer 
moves with a constant acceleration. Instead, the cold gas 
falls back quickly through the fingers. 

It is notable that th e char acteristic timescale for accel- 
erating the shell ^ ■y/2/o/a, with a ~ O.Olg is about 10^ 
years for the parameters under consideration. Since, the 
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Figure 11. Structures of the cold shell at times 2.3 X 10® yr (left 
panel) and 1.1 X 10^ yr (right panel). Initial optical depth across 
the cold shell is Ts = 1. Fingers and bubbles due to RTI are formed, 
which allows the cold dense gas to fall downwards instead of being 
pushed away by the background radiation flux. The vectors are 
the velocity field. 
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Figure 12. Top: History of the fraction of cold dense gas that is 
below the line -0.5«o- Time unit is 2.3 X 10® yr. Before the RTI 
develops, no cold dense gas is below the line —O.blg because the 
cold gas is pushed away by the radiation field. After RTI grows 
significantly, over 13% of cold gas falls back during ~ 2.8 X 10® 
yr. Bottom: Change of the maximum and minimum heights of 
the cold dense gas with time. During the time t < 3.6, the whole 
dense shell moves upwards. Once RTI grows significantly, some of 
the cold gas almost falls back to the bottom of the simulation box. 



timescale for appreciable growth of the RTI is a few x 10^ 
years when Ts ~ 1, there is sufficient time for the RTI 
to grow and disrupt the sheh. We have also performed 
runs with the radiation force is double the gravitational 
force, but keeping ~ 1. This increases the effective ac- 
celeration by a factor of ~ 100 and reduces the timescale 
for acceleration to ^ 10^ years. In this case the shell is 
accelerated efficiently and reaches the upper boundary 
of the domain before the RTI has time to grow appre- 
ciable. Hence we find that efficient acceleration of gas 
requires low optical depths or very super-Eddington ra- 
diation fluxes. 
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Figure 13. The same as Figure [TTI except that the initial optical 
depth across the cold shell is Ts = 0.01 in this case. The two 
interfaces are both stable in this case and the whole cold shell is 
pushed upwards. 



6. SUMMARY AND DISCUSSIONS 

Using our recently developed radiation hydrodynamic 
module for Athena, we have performed a set of numeri- 
cal simulations to see the effects of strong radiation field 
on the development of the RTI for a background state 
with pure scattering opacity. In many respects, the ra- 
diative RTI is rather similar to the purely hydrodynamic 
RTI. Instability is always present when a high density 
fluid overlies a lower density fluid and we find growth 
rates that are within an order of magnitude of the hy- 
drodynamic case, even when radiation pressure exceeds 
gas pressure by several orders of magnitude. Neverthe- 
less, we find that the presence of radiation modified the 
development of the RTI in several significant ways. 

For a gas pressure supported interface, the presence of 
a radiation field will generally reduce the growth rate of 
RTI, and growth rate decreases as the radiation pressure 
increases. This is because radiation acts as drag force 
when the material tries to move with respect to the radi- 
ation field. When radiation pressure is large (P ~ C), the 
development of small scale structure is also suppressed 
by the radiation held. 

We find similar effects for a radiation pressure sup- 
ported interface. In optically thick limit with an isotropic 
radiation field, the growth rate of the radiative RTI is 
similar to the non-radiative case. As we increase the ra- 
diation pressure and the opacity we flnd lower growth 
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rates, which we again attribute to radiation drag effects. 
If the optical depth across the domain is low, so that 
the radiation becomes anisotropic, we generally find the 
growth rate is reduced relative to the isotropic case for 
the same parameters. 

For a radiation supported interface with random per- 
turbations, we find that it is easier for the dense fingers 
to sink than the light bubbles to rise up in contrast to 
the hydrodynamic RTI, where rising and sinking is nearly 
symmetric. This effect is more significant when opacity 
is decreased and the diffusion time across the fingers is 
reduced. The mixing process is also suppressed in ra- 
diation RTI, especially when radiation pressure is much 
larger than gas pressure. 

The degree of radiation damping depends both on the 
opacity and the ratio between radiation pressure and gas 
pressure. When radiation pressure is large (P ^ C), 
all the modes will be affected by the radiation damp- 
ing effect significantly as long as photons and gas are not 
decoupled for this mode. Turbulence due to RTI associ- 
ated with secondary instabilities (e.g. Kelvin-Helmholtz) 
is significantly suppressed. When radiation pressure is 
just comparable to the gas pressure, only the small scale 
structures will be damped. 

We have used a background state with pure scatter- 
ing opacity because a background state with non-zero 
absorption opacity and an interface with a density dis- 
continuity cannot generally obey both mechanical and 
thermodynamic equilibrium. This configuration will be 
particularly problematic in optically thick regime and 
the background state will evolve quickly because of the 
thermalization process. Any results obtained for systems 
with absorption opacity (e.g. lJacauet fc Krumholal2011|) 
will only be valid when the growth time of RTI is suffi- 
ciently shorter than the thermalization time or any other 
short timescale associated with the interface evolution 
(e.g. the evolution of an ionization front in an HII re- 
gion). 

In order to assess the effectiveness of radiative driving 
of cold gas outflows we simulate the acceleration of a cold, 
dense shell, for which radiation forces exceed gravity. We 
find that when the optical depth in the cold shell is ~ 1 
and the radiation field is just slightly above Eddington, 
the RTI can develop quickly and prevent the accelera- 
tion of an appreciable amount of cold gas. However, if 
the radiation field is substantially super-Eddington, or 
if cold gas is very optically thin so that the growth rate 
of RTI is significantly reduced, the shell can be acceler- 
ated efficiently before the onset of the RTI. These simu- 
lations are too simplified to draw definite conclusions on 
the impact of RTI on radiation feedback in star forming 
environments, but suggest that the RTI is likely to be 
important in some cases. Those simulations also imply 
that if the growth time scale of RTI is smaller than the 
typical acceleration time scale, RTI is likely to change the 
structure of the system significantly. As we find that the 
growth rate of RTI is similar to normal RTI when radi- 
ation pressure is not significant (P ^ C) and the growth 



rate will be reduced when radiation pressure is very sig- 
nificant (P ^ C), our results give a general criterion on 
when RTI will be important when a specific astrophysical 
system is considered. Future work will benefit from the 
use of more realistic dust opacity and the use of non-grey 
(i.e. multiple frequency bins) radiative transfer to differ- 
entiate the impact of UV, mid-infrared, and far-infrared 
photons. 

One caveat of our simulations is that they are done 
in 2D. It is known that RTI leads to different satura- 
tion states in 3D compared with 2D (e.g. Marinak ct al] 
Il995f ). We plan to explore these effects on specific prob- 
lems in future work. However, the linear growth rate 
is expected to be unchanged in 3D. We expect that the 
qualitative effects we find for the radiation field on the 
growth of the RTI, such as the importance of radiation 
drag in very radiation dominated regimes, will also be 
unchanged in 3D. 

Other avenues for future research include studying the 
effects of absorption opacity on radiation RTI. If the ther- 
mal time scale is much longer than the growth time of 
RTI, we anticipate that absorption opacity will have mi- 
nor effects on the growth and saturation of RTI. When 
the thermal time scale is much shorter than or compa- 
rable to the growth time of RTI, the effective equation 
of state of the fluid will be modified. Furthermore, the 
background state will also change dramatically in this 
case, which can also affect the saturation of RTI. The 
effects of absorption opacity on the RTI in specific astro- 
physic environments (e.g. with dust opacity in ULIRGS) 
is a topic of ongoing research that will be presented in 
future work. It is also important to include magnetic 
fields and examine the effects of radiation on magnetic 
RTI. Our radiation MHD code can also be used to study 
RTI in more specific astrophysical systems, such as H II 
regions with massive star formation, photosphere of ac- 
cretion disks and supernova remnants. 
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APPENDIX 
LINEAR ANALYSIS 



We present a simplified linear analysis of the Rayleigh- Taylor Instability in the case of a single interface between 
two constant density media with radiation, but in the limit of negligible absorption opacity. The background for this 
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problem is described in Section |3] and specified by eqs. ([7]) - ([TU)) . We will consider the case with purely radiation 
support (a = 1) so the both density and gas pressure are constant within each half plane. The background is assumed 
to be static with Fr = FqZ and absorption opacity as = Kp. With these assumptions, hydrostatic equilibrium reduces 
to PFqk — g. We consider the 2D problem with a background that is uniform in x. Variables above and below the 
interface will be denoted by subscripts + and — , respectively. We focus on the case with a diagonal Eddington tensor: 
fxx = fzz = 1/3 a nd fxz = 0. Note that thi s assumption precludes radiation viscosity although still allows effects due 
to radiation drag (jMihalas fc MihalaslflQSl . 

We begin by linearizing eqs. (|4]) with perturbations of the form f{x, z, t) — A{z) exp [i{kx — ujt)]. We denote the x 
and z components of the velocity by u and w, respectively. We denote the (scalar) radiation pressure by Pr and the 
components of x and z components of F by F^ and F^ , respectively. Constant background quantities are denoted by 
a subscript zero, e.g. background density and vertical flux are denoted by Fq and po. We then have: 



iijj6p + p iku + 



dw 



-iojpu + ika^Sp H — FagSFx + . 



2 '^^P 

-lujpw + a — — 
dz 



SPr + ikSF, 



•asSF, + 
dSF, 



C 



-SF, 



— 5F,+ik5P, 

dSPr 



dz 
aJFrr 



dz 



agSFz + nFoSp - 



dz 
C 

Aa^Pr 

" 1 

C 

AasPr 

c 

AOsPr 



0, 


(Al) 


0, 


(A2) 


0, 


(A3) 


0, 


(A4) 


0, 


(A5) 


0. 


(A6) 



Since the radiation field does not exchange any thermal energy with the gas when aa — 0, gas energy equation simply 
reduces to a statement that gas pressure perturbation are adiabatic SP — a^Sp, with ~ jPo/po- Note that since 
both Pq and po are constant, a is a constant. 

We further simplify the problem by selectively dropping terms of order C^^. In particular, we drop the time 
derivative terms in eqs. (|A4p - (|A6p and the fourth term on the rhs of eq. (|A4|) . These terms are always small relative 
to other terms when C is large (i.e. the fiow is non-relativistic). We nevertheless retain the last terms on the rhs of 
eqs. (|A5P and (jA6P because they include Pr, which can be much larger than Fq in the limit of large optical depth. It 
is useful to first define a characteristic frequency 

_ iPrg _ 4:PrKP 

" ~ IhC ~ C ' 

noting that this frequency generally varies with height, due to the linear dependence on the background Pr on z. 
We can use equations (IA2p and (|A5|) to solve for u and SF^ , obtaining 



and 



SF, 
Fo 



_ k /gSPr 

LO \(TsFo 

k{ijj + 



5Pr 



va^k'^ Sp 



LU GsFq guj'^ p 

Inserting these into the above equations yields (after some algebra) a system of equations 

dX 

dz 

with 



(A7) 



(A8) 



AX 



and 



A = 



X^^ 







SPr 5Fz Sp ^ 

Fq Fq Po 



-i^ 



V 



n 9 

2 ^ 











(A9) 



-1 / 



Here we have replaced w by the Lagrangian displacement using the definition w = —ico^. These equations represent 
for ODEs, but with non-constant matrix coefficients i.e. the terms with v vary with z. In what follows we will ignore 
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this dependence and assume constant coefficients. Since Pr — Pq — 9Poz, v is approximately constants for lengtliscales 
I <C PPo/(ffPo), Hence, it is a reasonable approximation for the subset of simulations with P ^ 1. 

With the above assumption, solutions take the form Xj_ exp (A^^z) with four eigenvalues A!^. The eigenvalues for this 
system are given by 

+ iv) " 



A^ - /fc" A" - fc 



= 0. 



The first set of eigenvalues A^ = have eigenvectors with (5P = a^(5p = and 



Fo 



^SPr = %. 



Uj{u! + w) 



These correspond to incompressible modes. The other set of compressible eigenvalues have A ± ^Jk? — u]{ijJ + iv)/a^ 
and has eigenvectors with 



5P- 



^dPr = 



puj^{ijj + iu) 
gk?' + w^A ' 
p{gijj^\ — iuj^i/) 
gk^ + uj^X ' 
k'^ujiiu + iv) 



-Fo 



gk^ + Lu'^X 



We now apply standard boundary conditions, the first being that the perturbations vanish far from the interface. 
The second set of conditions is given by matching solutions across the interface. In addition to the continuity of 
w, the non-trivial matchin g conditions come from integrating the vertical conservation relations in eqs. As 
iJacauet fc Krumhold (|2011f) discuss, these are equivalent to the continuity of the Lagrangian perturbations across the 
interface, which for our background correspond to 

PAP^ = V5Pr - pg£,, 
AP = 6P, 
AF, = SF,. 

Defining perturbations above and below the interface proportional to 

X± exp (Tkz) + X± exp (tA±2:), 

where X± correspond to the compressible eigenvalues and are both defined to have positive real parts. These matching 
conditions define a matrix — with X'^ — {x+jX+jX-iX-)^ 



I 



A = 



where we have defined 9± 
yields 



" 1 Wi 



\\gk 



1 

P+. 

s+ 
gk^^ 



-1 


-1 



-1 



\ 



gk^ ± n?X± and n = —iuj. The dispersion relation is obtained by setting det A = 0, which 



(A+p_ + Ap+) 



gk 



+ =0 



The quantity in the first set of parentheses has a real part which is always greater than zero. Hence if we look for 
purely imaginary solutions with growth rates n = — iw, we obtain 



gk- 



P- 



P++ P- 



which is equivalent to eq. (jl4p . the standard result for the inviscid, incompressible hydrodynamic problem. 

The robustness of this form for the dispersion relation may seem surprising, but is straightforward to understand. 
First note that the solution corresponds to the case where = so the solution on either side of the interface 
is a purely proportional to the incompressible eigenvector. This result stems from the fact that only the purely 
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incompressible eigenvectors permit the continuity of both SFz and 6P across the interface. In the limit that Sp ~ 
and dropping the same terms as above, eqs. (jAip - (jA6p reduce to 

dw 

iku H — — = 0, 
dz 

-iujpu + ik¥SPr = 0, 

-iujpw + P—-^ = 0, 
az 

which resemble the hydrodynamics case, but with gas press ure perturbations rep l aced by radiation pressure. This 
should be contrasted with the optically thin case discussed bv lJacquet fc Krumholzl ()2011[ ). where there is no SPr since 
the radiation field is decoupled from local conditions near the interface. In that case, since gravity is exactly balanced 
by radiation forces the effective gravity is zero and the interface is nominally stable. 

Finally, we note that this result seems to imply that radiation drag have no effect on the linear growth rate. However, 
when ly is large {v ^ Vgk), radiation drag forces make the linear regime of growth be rather limited. For large values 
of v, non-linear effects can become important, even for very modest velocity amplitudes. Indeed our numerical 
simulations indicate that non-linear effects become important at lower amplitudes for parameters corresponding to 
large v, as discussed in section 1411 
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